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TECHNICAL MEMORANDUM X- 64867 


ENHANCING SIMULATION EFFICIENCY WITH ANALYTICAL TOOLS* 

INTRODUCTION 


The relatively recent development of computers — particularly digital 
machines — has strongly influenced engineering design. Modern engineers and 
designers place heavy reliance on digital, analogue, and hybrid computers to 
simulate models of systems being designed. This is particularly true in the 
area of dynamics and control, where computers are used extensively to simulate 
the dynamic behavior of such models. In addition, digital computers often are 
used to solve the equations of motion describing the model of the dynamics of a 
particular system for the associated eigenvalues. 

An older (but continually updated) approach to system design is to per- 
form mathematical analysis of the representative equations of motion. This 
usually leads to closed form solutions of these equations. 

Recently NASA has placed particular emphasis on the development of low 
cost aerospace systems. The area of computer simulation has been identified 
as one of the many areas where cost reduction should be investigated. The 
purpose of this paper is to indicate some means of combining both computer 
simulation and analytical techniques in order to mutually enhance their efficiency 
as design tools and to motivate those involved in engineering design to consider 
using such combinations. While the idea is not new, heavy reliance on computers 
often seems to overshadow the potential utility of analytical tools. Although the 
example used in this paper is drawn from the area of dynamics and control, the 
principles espoused are applicable to other fields. 


APPROACH 


Computers provide a relatively rapid means of solving the equations of 
motion of models of physical systems. They can handle large order systems of 


* Portions of this paper were presented at the Eighth Annual Allerton Confer- 
ence on Circuit and System Theory at Montieello, 111., Oct. 7-9, 1970, and at 
the First Annual Research and Technology Review, George C. Marshall Space 
Flight Center, Alabama, Feb. 22-23, 1973. 



equations of a nature and dimension that historically have troubled analysts. 
These outstanding and well known traits often lead to exclusive reliance on 
computer simulation techniques by contemporary engineers. This exclusive 
reliance can lead to technical pitfalls in the form of erroneous results because 
of programing errors or to improper attention being devoted to implicit com- 
puter characteristics such as accuracy (round off errors) and clock time. 
Furthermore, suitable numerical values for parameters used in the simulation 
often are determined by ’’cut- and- try” methods, which can be time consuming 
and inefficient. 

Mathematical analysis, on the other hand, usually requires relatively 
low order systems with the attendant simplifications leading from the original, 
more complex (and more accurate) system model. Further, the user of 
present analytical techniques often encounters difficulty in dealing with non- 
linear systems or time varying parameters. However, when a relatively 
simple system model can be obtained, analytical techniques exist that can 
predict the dynamic behavior of that model, often as a closed form solution. 


From the foregoing, it is seen that analysis techniques and computer 
simulations might be used to mutually enhance their strengths and decrease 
their weaknesses. Obvious examples would be to use analytical techniques to 
predict the dynamic behavior of the computer simulation and thereby help ”de- 


-l-'U ^ 
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Another use is to apply analy Ileal techniques to determine 


numerical values for system parameters that are needed for the computer 
simulation, thereby reducing or eliminating ”cut-and-try” efforts. It is 
apparent that the use of simplified models might lead to erroneous results. 

This danger is reduced by comparing the dynamics of the simplified model 
used for analysis with the dynamics of the computer simulation. Disagreement 
indicates an error in the choice of the simplified model, a mistake in the 
associated analysis, or a programing error. While such disagreements are 
inconvenient to resolve, once an agreement is achieved the simulation results 
may be used with a high degree of confidence. The major portion of the 
remainder of this paper will be devoted to showing how an analysis technique 
may be used efficiently to select numerical values for a large scale computer 
simulation. This was actually performed in the design and development of the 
Sky lab control system. 


SKYLAB EXAMPLE 


In late 1970 a hybrid simulation of the Sky lab control system, simulating 
the rigid body dynamics, was under development for use in the design of the 
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control system being developed for NASA’s first manned orbital space station, 
Skylab (Fig. l)[ lj. It was found that the system response was unstable unless 
the sampling period was decreased to an unacceptably small value in terms of 
computer utilization. Hence, a simple digital filter was introduced into the 
system. The effect of the numerical values chosen for the basic control sys- 
tem parameters and filter parameters on the system dynamics, and the sampl- 
ing period, was analyzed and is described below'. 



Figure 1. Skylab. 


Definition of Problem 

The Skylab control system is shown in block diagram form in Figure 2. 
If it is assumed that the reaction moment (M ) is kept equal to the commanded 

moment ( M^) and only rotation about a single axis is considered, the block 

diagram may be simplified as shown in Figure 3. If at first no digital filter is 
present (K 2 = 0), the problem is to choose values for control gains K 0 and K lf 
and sampling period (t) such that desirable transient characteristics accrue. 
However, as will be shown subsequently, a simple digital filter (indicated on 
Fig. X) must be added to the system to achieve acceptable sampling periods. 
Now, the problem becomes one of choosing values for the control gains K 0 and 
K t , and the filter element K 2 so that an acceptable value of T will provide 
desirable transient characteristics. 




3 


1 


VEHICLE 

ATTITUDE 

KEP.«<$ r > 



gi&sbal 

RATE 

ERftOH 


Figure 2. Block diagram of Sky lab control system. 


Description of Parameter Plane Method 

The parameter plane technique for analysis and synthesis of linear and 
nonlinear control systems is amply described in Siljak f s recent monograph on 
the subject [2]. Reference 3 describes the application of the technique to the 
analysis and synthesis of linear sampled-data control systems. 



Figure 3. Simplified control system block diagram. 

Once the system characteristic equation has been obtained, the param^ 
eter plane method enables the designer to evaluate graphically the roots of the 
equation. Hence, he may design the control system in terms of the chosen 
performance criteria; e.g\, absolute stability, damping ratio, and settling 
time. He is able to see the effect on the characteristic equation roots of 
changing two adjustable parameters. Siljak further simplified the design 
pioceduie by introducing Chebyshev functions into the equations, thereby 
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putting them in a form particularly amenable to their solution by a digital 
computer. The method has been extended to portray the effect of varying the 
sampling period [4]. Thus one can see the effect of the choice of values 
assigned to the sampling period on absolute and relative stability. Also, the 
recursive formulas shown herein are simpler in form than the Chebyshev 
functions [3]. Although the brevity of this paper precludes an explicit demon- 
stration, the resulting formulation is deliberately cast in a form that makes 
it particularly amenable to solution by a digital computer or a desk calculator,, 
again emphasizing the interplay between analysis and computing machines. 
Because a calculator is usually readily available to the design engineer, the 
example portrayed here was solved using a Hewlett-Packard desk calculator. 

The technique requires that the control system be described by a 
characteristic equation which is transformed into the z-domain. Two adjust- 
able parameters ( k 0 , k t ) are selected, and the characteristic equation (CE) 
is recast in terms of them; i.e., 
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(1) 

y. « a.k 0 + b.k. + c. , 

J J J 1 J 


(2) 

Ts id 

z = € = re , 


(3) 

T 

r = r ^’V T) = e 11 

y 

(4) 


and 

/3 = T) = cos 9 = cos (cu^Tn/i- £ 2 ) , ( 5) 

where £,gj , and T represent the damping ratio, natural frequency, and 
sampling period, respectively. To transform the characteristic equation from 
a differential equation into an algebraic equation, z J may be defined as 
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z J = X. + iY. sll - (3 2 , (6) 

J J 


where X 0 = 1, Y 0 = 0, X t = r/3, and Y t = r. The following recursion formulas 

are found to enable one to solve for other values of X. and Y.: 

J ] 


X. - 2r {3X. + r 2 X. , = 0 
J+l J J- 1 


and ( 7) 


Y. - 2r£Y + r 2 Y. = 0 

J J-l 

If equation (6) is substituted into equation ( l) and the real and imaginary 
parts of the resulting equation are separated, two simultaneous algebraic 
equations are obtained* They may be solved for the adjustable gains k 0 and 
k t : 


kn ~ 


BiC 2 - B 2 Ci 


k, = 


A 2 Ci - A 1 C 2 


( 8 ) 


J - AjB 2 - A 2 Bji 


( 9 ) 


A i = V. a.X. 
.11 
3 = 0 J J 


B t = y b.x. 

• ii 

j=0 J J 


c t = f. c.x. 
j-0 j J 


n 


A 2 = 2 a i Y i 


j = 0 J J 


Bj = Y b.Y. 
t - J 1 1 
j=o J J 


c 2 = f. c.Y. 

. 1 i 

3=0 J J 


( 10 ) 
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Recalling that and Y, are functions of and T and observing that a^, 

b , and c. may be functions of T, we see that contours for specified values of 
J J 

?> may be plotted as functions of oj and T in the k 0 - k t parameter plane. The 


locations of these contours in the parameter plane represent values of the 
complex conjugate roots of the characteristic equation. If one chooses 
numerical values for and values for parameters and and the 

sampling period T are established. If the coefficients a., b., and c. contain 

J J J 

exponential terms with T appearing in the exponents, then each exponent must 
be replaced by its power series and truncated according to the accuracy that 
is desired. 


The real root locations corresponding to values of z when 9 equals 
0 degrees and 180 degrees may also be plotted on the parameter plane by set- 
ting z equal to a positive or negative real constant, substituting that value 
into equation ( 1) , and solving for as a function of k 0 . Each resulting contour 
corresponds to a location of a real root in the z-domain. When z = +a, 


2 V = 0 ’ ( 1J ) 

3=o 1 


and when z = -a, 


a/2 .. 

Z ? 2j a 

3 = 0 J 


( n-2) / 2 


j = 0 


Y (2j + l) a 


(2j + l) 


, n even , 


( 12 a) 


( n_ ^ /2 

j-0 


2j (2j+l) 

y a - y, a 
2j (2j+l) 


, n odd , 


(12b) 


where a is a positive real number. 

For a linear sampled data control system to be stable, it is necessary 
that all roots of the characteristic equation lie within the unit circle on the 
z-plane. It is assumed that the system being designed possesses low-pass 
filter characteristics so that only the primary strip ( corresponding to 0 < 0 < 7 r 
in the z-plane) need be considered. The stable region is bounded by the semi- 
circle defined by the upper half of the unit circle. This region may be defined 
by mapping three contours from the z-plane onto the parameter plane. The 
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first contour' is defined by setting £= 0 in equations (8). The resulting complex 
root stability boundary is seen to be a function of co and T. The second and 

third contours correspond to real root locations at z = +1 and z = -1 and are 

found by setting "a" equal to unity in equations (ll) and (12), respectively. 

The stable region is determined by applying a shading criterion or using a test 

point [3], If the Jacobian, defined by equation (9), is greater than zero, then 

the stable region ( if it exists) lies to the left of the £ = 0 contour as u> T 

- n 

increases; the left side of the line is double crosshatched to indicate a boundary 
associated with double, or complex conjugate, roots, (if the Jacobian is less 
than zero, the stable region lies to the right.) Single crosshatching is used 
on the two contours associated with the real roots. The side of the contour 
on which to place the erosshatching is determined by the requirement that 
crosshatching be continuous, or on the same side, of the contours as the 
intersections corresponding to z = +1 and z = -1 are approached along either 
the complex root or real root stability boundary. 

Once the stable region, if it exists, has been delineated in the param- 
eter plane, the transient characteristics of the system can be determined in 
terms of the locations of the roots of the characteristic equation. For the 
complex conjugate roots, these locations are defined in terms of damping 
ratio and natural frequency. Contours of constant L are determined as func- 
tions of o> and T from equation ( 8) and plotted on the parameter plane. 

Similarly, values of "a" corresponding to real root locations may be sub- 
stituted into equations (ll) and (12) to plot contours corresponding to these 
locations. Now the effect of any chosen design point on the parameter plane 

may be found in terms of C, go , a, and T. 

n 

The analytical technique developed permits the designer to observe the 
effect of simultaneously changing three control parameters and the sampling 
period. Most existing conventional techniques permit the observation of the 
effect of changing only one control parameter and do not show the effect of 
various sampling periods. 

Sometimes it is specified that the settling time of the system be less 
than a prescribed value. This corresponds to requiring that the real part of 
the roots of the characteristic equation be less than a prescribed negative 
real constant. A boundary corresponding to this requirement can be drawn 
on the parameter plane by mapping a circle of constant radius (for a chosen 
value of c o t, and T) from the z-plane onto the parameter plane [3]. Relations 

exist for estimating the maximum overshoot and peak time of transient 
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response when it is valid to assume a second order system [5]. However, 
a simple estimate can sometimes be made by merely looking' at t^ie difference 
equation representing the system response, estimating when the overshoot 
will occur, and plotting a corresponding line on the parameter plane. This 
procedure is brought out in the example. Steady state response may be found 
from the open loop transfer function and the assumed forcing functions in the 
conventional manner [ 5] . 


Application to Skylab Design Problems 

If at first no digital filter is present, the problem is to choose values 
for control gains k 0 and k t and sampling period T such that desirable transient 
characteristics accrue. The open loop transfer may be written in the z- 
domain: 


G( Z ) = * j (* + Kl ) (M~ )( e " T d s )(ir4 ( 13) 

If the computation delay T^ is made equal to the sampling period T, then equa- 
tion ( 13) may be rewritten as 

ni„\ - ( k o + k t ) z + ( k 0 - lq) 

Z (Z-1)2 ’ < 14 >. 


where 


ko 


K 0 T 2 

21 


and 


k i 


K t T 

I 


( 15) 
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The corresponding characteristic equation is equation (l), where n = 3. Table 

1 defines the numerical values of the coefficients a., b.. and c for this case 

J J J 


if K 2 is set equal to zero. 


TABLE 1. COEFFICIENTS OF CHARACTERISTIC EQUATION 


j 

a. 

J 

b. 

J 

c. 

J 

0 

1 

-1 

k 2 

1 

1 

1 

1 - 2K 2 

2 

0 

0 

K 2 - 2 

CO 

0 

0 

1 


Using Table 1 values and setting £ = 0 and a = 1 for equations (4), (5), (7), and 

(8) through (l2) provides J<0 and the three stability boundaries shown in 

Figure 4 (again setting K 2 = 0). The number of stable roots in each region is 

indicated in parentheses. It is seen that the stable region lies in the first 

quadrant. Contours corresponding to values of Z, >0 may now be plotted as 

functions of w T. Several are shown on Figure 5 with equal values of to T 
11 n 


connected by dashed lines. 


K, 



Figure 4. Stability contours with filter. 
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Figure 5. Relative stabilit}' contours (original system, K 2 = 0). 

If a value of £ = 0. 7 is desired, o^T must be less than 0. 8. The predicted 

system natural frequency is 0. 17T, yielding a maximum possible sampling time 
of approximately 2.4 seconds. To obtain acceptable transient characteristics, 
the actual value of T used would have to be considerably less than 2.4 seconds., 
With the goal of utilizing less computer time, it is desired to add a simple 
digital filter to increase the sampling period. A filter of the form, 


F(z) 


z 

z+ K 2 


(16) 


is introduced (Fig. 3). The qualitative effect of the filter may be seen by 
referring to Figure 4. For positive values of K 2 , the z = -1 stability boundary 
is displaced upward. The effect is to move the £ = 0 curve so that it bulges 
more to the right and upward, enclosing a larger region and increasing the 
maximum value of w T for stability. An example is shown for K 2 = 1 (Fig. 6). 
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Figure b. Relative stability contours ( K 2 = 1). 

Howcvei , if the value of K 2 is set too large, the z — — 1 boundary will cut off 
part ( or all) of the stability boundary, as indicated in Figure 7. As K-> assumes 
negative values, the z = - 1 boundary is displaced downward and the £> = 0 curve 
moves to the left until, for values of K 2 < - 1, no stability region remains. The 
steady state error is zero for step displacement and velocity inputs. For 
acceleration inputs, it is a constant value which is proportional to Kr> and 
inversely proportional to k 0 . 

First a value of K 2 should be chosen. It is seen that as K 2 increases 
positively, larger values of co^T will lie in the stable region for corresponding 

values of Although the purpose of the filter is to increase the value of T, 

K 2 mLls t not be increased to the value where the stable region disappears or to 
a value greatei than the computer clock time. The clock time is 2. 5 milli- 
seconds, and the simulation runs are performed at 100 times real time, so 
the sampling period T must not exceed 2. 5 seconds (go T< 0. 79). Looking at 

the difference equation associated with the system response 0, one sees that 
in geneial large values of k 0 and kj will result in large overshoots of the 
response: 
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0(nT) = (k 0 + k,) 0 r [(n - 2)T] + ( k 0 - k,) ^ [ (n - 3)T] 

-(K 2 -2) <f> [(n- 1)T) -(ko+ k,+ 1 - 2Kj) <p [(n - 2)T) 

-( ko - k[ + K 2 ) <p [ ( n - 3) T) . ( 17) 

As a design compromise, a value of K 2 = 1 was chosen. If a value of 

0. 707 is used for £, the maximum value of to T that will remain within the 

n 

stable region is 2. 1, which is greater than the clock-time constraint. Choos- 
ing values for co^T on a particular ^.-contour establishes the location of the 

real root. 


13 



If both the difference equation ( 17) and the parameter plane contours are 

examined simultaneously, one usually can quickly deduce the sampling period 

in which the maximum overshoot will occur. For example, if it is assumed 

that the unit step response will have a maximum overshoot in the vicinity of 

the third sampling instant, contours corresponding to selected magnitudes of 

the response at the third sampling instant may be mapped onto the parameter 

plane, using equation (17). Contours corresponding to unit step response 

values at the third sampling instant of 1.2, 1.4, and 1.6 are indicated on 

Figure 6. If it is desired that the overshoot not exceed 60 percent, then it is 

estimated that 1.5 <co T< 0. 7. Because of the clock-time constraint, T is 

a n 

set at the lower value. The resulting step input response is shown on Figure 8. 


0 



Figure 8. Unit step input response. 
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OTHER EXAMPLES 


Other examples are available which further amplify the prudence of not 
accepting computer results without checking them by analysis, Schiehlen [6] 
found that the results reported in Reference 7 are incorrect ' — probably due 
to round-off errors in the digital computation. Similar pitfalls exist with 
analogue computation: the difficulty in obtaining the well known stability bound- 
aries of the Mathiew equation is usually a frustrating exercise in futility [6]. 

CONCLUSION 


The main body of this paper describes an analytical technique that was 
used to predict desirable numerical values for simulation parameters. Numer- 
ical values were rapidly found and used successfully in the hybrid simulation. 
This combination of an analytical technique with the simulation technique 
obviated the need for "cut-and-try” methods to choose the numerical values, 
thereby saving both time and computer utilization. 
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